function [obj,grad,var] =  LikNPL(theta,F,States,StackedAction,StackedSingleState,IndFieldID,pproc,p,opt,do_var)
    
Ntheta = length(theta);
Regions = unique(States.RegionsForExoStates);
%NumR = length(Regions);
StackedSize = size(StackedAction,1);
Nstates = size(States.AllStatesMap,1);

TU = opt.Years;
T = length(TU);
StackedRegionId = kron(ones(T,1),States.RegionsForExoStates);
dlogPitdTheta = zeros(StackedSize,Ntheta);
logPit        = zeros(StackedSize,1);

StackedA0 = (StackedAction == 0);
StackedA1 = (StackedAction == 1);
StackedA3 = (StackedAction == -1);

for r = Regions'
    reg = r;
    [U{r},dUdTheta{r}] = Payoff(States.AllStatesMap,States.MidPointsFixedStates,States.FixedStatesMap,States.SRDMidPoints,theta,pproc.y{reg});
    % Renewal choice we will work with is choice 1 (Start sugar or Replant)
    
    %load([opt.ScratchFolder 'F' num2str(reg) '.mat'],'F');
    
    % Compute Log(P): rows -> states, column -> Action.
    for a=1:3
%        ExpA(:,a) = exp(U{r}(:,a) - opt.rho*F{r}{a}*log(p{r}(:,2)));
        %ExpA(:,a) = exp(U{r}(:,a) - opt.rho*F{r}{a}*log(p{r}(:,2)));
         ExpA(:,a) = exp(U{r}(:,a) + opt.rho*F{r}{a}*( U{r}(:,2) - log(p{r}(:,2))) );
    end
    
    %clear F
    
    DenA = sum(ExpA,2);
    LogP{r} = log( bsxfun(@ldivide,DenA,ExpA) );

    % Compute dLog(P)dTheta
    % First there is a term common to all choices for dLog(P)dTheta
    for a = 1:3
%        Aux1{a} = bsxfun(@times,ExpA(:,a),dUdTheta{r}{a});
         Aux3{a} = opt.rho*F{r}{a}*dUdTheta{r}{2};
         Aux1{a} = bsxfun(@times,ExpA(:,a),dUdTheta{r}{a} + Aux3{a});
    end
    Aux2 = Aux1{1} + Aux1{2} + Aux1{3};
    Aux2 = bsxfun(@ldivide,DenA,Aux2);
    for a=1:3
%         dLogPdTheta{r}{a} = dUdTheta{r}{a} - Aux2;
         dLogPdTheta{r}{a} = dUdTheta{r}{a} + Aux3{a} - Aux2;
    end

    % Now compute the matrix dPitdTheta. Element (n,m) is the partial
    % derivative of probability of chosing the choice made at observation n at
    % the state of observation n wrt theta m.

    Action0r = (StackedA0 & StackedRegionId == reg & StackedSingleState > 0);
    Action1r = (StackedA1 & StackedRegionId == reg & StackedSingleState > 0);
    Action3r = (StackedA3 & StackedRegionId == reg & StackedSingleState > 0);
    sA0r = StackedSingleState(Action0r);
    sA1r = StackedSingleState(Action1r);
    sA3r = StackedSingleState(Action3r);

    dlogPitdTheta(Action0r,:) = dLogPdTheta{r}{1}(sA0r,:);
    dlogPitdTheta(Action1r,:) = dLogPdTheta{r}{2}(sA1r,:);
    dlogPitdTheta(Action3r,:) = dLogPdTheta{r}{3}(sA3r,:);

    % Save log P for objective function:
    logPit(Action0r) = LogP{r}(sA0r,1);
    logPit(Action1r) = LogP{r}(sA1r,2);
    logPit(Action3r) = LogP{r}(sA3r,3);
end

% For the next step we need to sum the above by observation across time:
% SumTdlogPitdTheta(n,m) will be the sum of partial der wrt theta_m for field n across
% time.

NumFields = StackedSize/T;
SumTdlogPitdTheta = zeros(NumFields,Ntheta); % allocate space.

for j=1:length(theta)
    SumTdlogPitdTheta(:,j) = accumarray(IndFieldID,dlogPitdTheta(:,j));
end

obj  = - sum(logPit);
grad = - sum(SumTdlogPitdTheta,1)';

if do_var == 1
    %NAc  = sum(~isnan(DSe.Action));
    var  = ((1/NumFields)*(SumTdlogPitdTheta'*SumTdlogPitdTheta))\eye(size(theta,1));
    var = var/NumFields;
else
    var = nan;
end
    
    

    
end